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Abstract. In this paper we study some aspects of search for an immobile target by a 
swarm of N non-communicating, randomly moving searchers (numbered by the index 
k, k = 1, 2, . . . , Af), which all start their random motion simultaneously at the same 
point in space. For each realization of the search process, we record the unordered set 
of time moments {t^}, where Tk is the time of the first passage of the k-th searcher to 
the location of the target. Clearly, t^'s are independent, identically distributed random 
variables with the same distribution function ^(t). We evaluate then the distribution 
P{uj) of the random variable w ^ ti/t, where r = A^^^ X^fcLi ''"fc the ensemble- 
averaged realization-dependent first passage time. We show that P{uj) exhibits quite 
a non-trivial and sometimes a counterintuitive behaviour. We demonstrate that in 
some well-studied cases {e.g., Brownian motion in finite d-dimensional domains) the 
mean first passage time is not a robust measure of the search efficiency, despite the 
fact that 4'(t) has moments of arbitrary order. This implies, in particular, that even 
in this simplest case (not saying about complex systems and/or anomalous diffusion) 
first passage data extracted from a single particle tracking should be regarded with an 
appropriate caution because of the significant sample-to-samplc fiuctuations. 
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Introduction 

Search processes are ubiquitous in Nature: In order to survive, predators have to hunt 
the prey and the prey have to forage [H [21 13] . In order to convert into required reaction 
products, the reactants involved in chemical or biochemical reactions have first to find 
each other [1]. In many biophysical processes ligands search for binding sites [HI Ej, 
proteins seek the target sequences on DNA's [7], etc. Human beings look for a better 
job, partners, shelter, files in databases. Even an attempt to unlock a pin-protected 
device can be considered as a search in the space of all possible passwords [8] . 

Search for a desired target depends generally on a variety of different conditions 
and may take place in different environments: targets may be sparse, hidden, difficult to 
detect even when found. The targets may be immobile or mobile, try to avoid searchers 
or to evade from the searched area P, [TD] . They may have no or may have their own 
life-time and vanish before they are detected. 

Searchers may be immobile, as it happens, e.g., in visual search [H], may move 
freely or interact with the environment. Their motion may be hindered under conditions 
of molecular crowding [12], e.g., in cell's cytoplasm or in dynamical backgrounds formed 
by other randomly moving particles, or facilitated due to interactions with the molecular 
motors. The searchers may search "blindly" detecting the target only upon an encounter 
with it, or "smell" (or "see") the target somehow at long distances correcting their 
motion [9l [13]. They may have no memory of previously visited area or adapt their 
strategy "on-line" repelling themselves from their footprints on the searched substrate. 
Finally, the searchers may act individually or in swarms [HI [151 [IB] . 

In general, for each specific situation different search strategies may be realized, 
and the question of efficient ones has motivated a great deal of work within the last 
years. While earlier works have considered deterministic search algorithms (see, e.g., 
Refs. [H [21 [17] and references therein) specific to such human activities as, say, search for 
natural resources or rescue operations, more recent studies focused on random search 
strategies. It was realized that the strategies based on Levy flights or walks [H [THl [l9] , 
in which a searcher performs excursions whose lengths are random variables with heavy- 
tailed distributions, in some aspects are more advantageous than a search based on a 
conventional Brownian motion, or on random walks which step on nearest-neighbors 
only. Naturally, in this case the large-scale dynamics of searchers is super diffusive. 

Following the observation of trajectories of foraging animals in which active local 
search phases randomly alternate with relocation phases (see, e.g., Refs. [20]) another 
type of random search - an intermittent search - has been proposed. In this algorithm 
the search process is characterized by two distinct types of motion - ballistic relocation 
stage when the searcher is non-receptive to the target and a relatively slow phase with a 
random Brownian- type motion when the target may be detected [211 [221 [221 [21]. Much 
effort has been invested recently in understanding different optimization schemes for 
such a random search. In particular, one looked for the conditions allowing to minimize 
the mean first passage time for the process which is unlimited in time [21], or seeked to 
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enhance the chances of successful detection by minimizing the non-detection probabihty 
for the search process constrained to happen within a finite time interval [221 l23l 123] . 
Note that for such a search the large-scale dynamics is diffusive, albeit intermittent. 

Finally, a combination of a Levy-based and intermittent search has been proposed in 
Ref . [25] , in which the length of the relocation stage was taken as a random variable with 
a heavy-tailed distribution. It was shown that such a combined strategy is advantageous 
in the critical case of rare targets. 

In this paper we discuss some aspects of a blind search by a swarm of independent, 
non- communicating searchers. We consider a situation, as depicted in Fig. [H in which 
searchers (numbered by the index k, k = 1,2, . . . , N) occupy initially the same 
position in space, at some distance Xq apart of an immobile target, start their random 
motion simultaneously and arrive for the first time to the location of the target at times 
Tfc, respectively. Note that r^s are not ordered. Clearly, in such a situation the first 
passage times r^'s are independent, identically distributed random variables with the 
same distribution \E'(r). 

We focus here on the random variable 

1\ T 

where r is the averaged, over the ensemble of searchers, realization-dependent first 
passage time, 

1 ^ 

k=l 

Hence, the random variable u probes the first passage time of a given searcher relative 
to the ensemble-averaged first passage time for independent searchers. The scaling 
factor 1/A^ in Eq. (IT|) is introduced here for convenience, so that regardless of the value of 
A^, the random variable u has a support on [0, 1]. We note parenthetically that random 
variables such as in Eq. ([T]) were previously studied in Refs. [26l |27l [28], [29], [30], [31] within 
a different context. 

Our goal is to calculate the distribution function 

P{u) = (6(u-l^]), (3) 



A^ 



r 



where is the delta-function and the angle brackets denote averaging over different 
realizations of independent, identically distributed random variables r^. 

For arbitrary normalized "^{t), the distribution P{u}) is normalized, possesses all 
moments and the first moment (u) = udu P{uj) = 1/N. Note that P{u}) can be seen 
as a measure of the robustness of a given search algorithm. Clearly, if P{u}) appears to 
be sharply peaked at a; = 1/A^, this would signify that the underlying search algorithm 
is quite robust. Otherwise, if the distribution appears to be broad, or even to have 
a multi-modal shape, this would imply that the performance of such an algorithm is 
rather poor and sample-to-sample fluctuations matter. From yet another conceptual 
perspective, one can say that P{oj) probes the validity of the mean first passage time 
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Target 



Figure 1. A sketch of trajectories of four searcliers starting at the same point (filled 
black circle) at a distance xq from the target and reaching the target for the first time 
at different time moments. 



as the proper measure of a search process efficiency. We proceed to show that in many 
situations P{u}) has a rather comphcated structure of which the mean behavior is not 
representative. We will show that in different situations P{uj) may have completely 
different shapes (modality) and also change the shape (say, from a unimodal bell-shaped 
form to a bimodal f/- or an M-shaped one) when some of the parameters are slightly 
modified. A similar phenomenon of a shape reversal has been previously observed for 
different mathematical objects in Refs. [27 t |30 | [3T|. 

We turn next to the key feature: the distribution \E'(r) of r^'s, which encodes 
all the information on the specific properties of the searchers' random motion and on 
their initial location relative to the target. Without a significant lack of generality, we 
suppose that for a search in infinitely large systems, in which all searchers perform a 
blind search and detect the target with probability 1 upon a ffist encounter with it, 
\E'(r) can be defined as 

*(^) - W) (-7) 7^- W 

where r(-) is the gamma function, a is a characteristic parameter which sets the cutoff of 
the distribution at small values of r, and /i > is the so-called the persistence exponent 
|32j . Note that depending on the dimensionality of space and type of random motion, 
one can encounter completely different values of fi; fi can beO</i<l,/i = lor/i>l. 
The distribution in Eq. is normalized but does not possess already a ffist moment, 
(i.e., the mean ffist passage time), for < /i < 1. 
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The distribution in Eq. (jl]) is exact for a Brownian motion (BM) in semi-infinite 
one-dimensional (ID) systems, in which case /i = 1/2 and a = x^/AD, D being the 
diffusion coefficient [33]. In fact, the case /i = 1/2 appears to be very representative. 
According to the theorem due to Sparre Andersen [31], in ID systems for any discrete- 
time random walk with each step length chosen from a continuous, symmetric but 
otherwise arbitrary distribution the first passage time distribution "^{n) decays with 
a number of steps n as n~^^'^. For continuous-time r Markov processes, an analogous 
result is ^I/(t) ~ r~^/^. This universality, of course, is broken when subordination effects 
(i.e., long tailed waiting time distributions in a CTRW sense) are present. 

Although not exact, the distribution in Eq. (jlj) is physically quite plausible for other 
types of random motion in infinite ID systems provided that an appropriate choice of 
the exponent fi is made [33]: In particular, Eq. (j4]) with /i = 1 — 1/aisa reasonable 
approximation for the first-arrival probability density for Levy flights with Levy index 
a, 1 < a < 2 [35] and with fi = 1 — H [36] - for the first-passage-time distribution for 
fractional BM with Hurst index H, < H < 1. One may also claim that fi = 1 — df/dw 
for df < dyj (compact exploration [37]), where df is a non-integer spatial dimension 
and dy^ is the fractal dimension of random motion trajectories, since the first passage 
distribution is defined as the time derivative of the survival probability of an immobile 
target or of the normalized current through the surface of the target [33]. Hence, its 
long-time tail is the same as the long-time tail of the time derivative of the Smoluchowski 
constant [38| 139]. 

The case /x = can be encountered in border-line situations of compact exploration, 
when the dimension of space df equals the fractal dimension d^ of the random motion 
trajectories. In particular, this situation is realized for standard BM [d^ = 2) in two- 
dimensional (2D) space. Here, one finds that the long-time tail of the distribution \E'(r) 
follows (see, e.g., Ref.[33]) 

*(^) ~ -r^T-y 

rln (r) 

This situation will be considered in more detail in what follows. 

Finally, one finds fi = d/2 — 1, (or, more generally, fi = df/d^ — 1), for Brownian 
(anomalous) motion in d > 2 [df > d^) dimensional systems. Note, however, that 
here the expression in Eq. (jl]) is not normalized and thus can not be considered as a 
distribution since only a finite fraction of trajectories will visit the target within an 
infinite time so that the searchers will have a finite probability to escape to infinity 
(non compact exploration [37]). One can, however, normalize ^I/(r) in Eq. (jl]) by hand, 
treating it as the conditional probability distribution of the first passage times for such 
trajectories which visit the target within an infinite time. 

Next, we will consider a normalized, exponentially-truncated version of the 
distribution in Eq. (jl]): 

= -J^^^^ exp (-^) 4- exp (-^ , (6) 
2K^{2y^) V rJ r^+^ ^ V 67 ' 
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where K^{-) is the modified Bessel function. In contrast to \I/(r) in Eq. (jl]), the 
distribution in Eq. (|6]) possesses moments of arbitrary order. Note that the latter point 
is crucial and, according to a common belief, the first moment of the distribution, i.e., 
the mean first passage time can be regarded as a robust measure of the search process 
efficiency. We set out to show that in many situations this is not the case due to the 
significant sample-to-sample fiuctuations. 

The distribution in Eq. ([6]) is exact for a BM in semi-infinite ID systems in presence 
of a constant bias pointing towards the target (see, e.g., Ref.[33]). In this case one has 
a = Xq/AD, where Xq is the starting point and b = AD/v"^, v being the drift velocity. 

One may argue, as well, that the expression in Eq. ([6]) is an appropriate 
approximation for the first-passage time distribution for search in finite systems or 
search assisted by smell. In general, of course, for random motion in finite systems the 
distribution ^'(r) will be represented as a series of exponentials; the form in Eq. ([6]) 
is thus tantamount to a heuristic approximation of this series in which one takes a 
behavior specific to an infinite system and truncates it by an exponential function with 
the characteristic decay time equal to the largest relaxation time. Depending on a 
particular situation, the cut-off parameter b will be either proportional to the volume 
of the system (for BM in d > 2 systems), to where L is the length of the interval for 
BM in ID, to S'ln(S') for BM in finite 2D systems of area S, or inversely proportional to 
the strength of the bias for biased diffusion. Some of these situations will be discussed 
in detail in the sequel. 

We finally remark that the distribution in Eq. ([6]) appears in many other physical 
problems. To name but a few we mention the distribution of times between action 
potentials (or the ISI distribution) in the integrate-and-fire model of neuron dynamics 
[lO], the distribution of the stopping distances for the sliding motion of a solid block on 
an inclined heterogeneous plane [H], the avalanche life-time distribution in the mean- 
field version of the Bak-Sneppen model [12], the distribution of the probability current 
in finite disordered one-dimensional samples |13] or the distribution of the number of 
times that a particle diffusing in a sphere hits its boundary during some time interval 
[6]. Thus our subsequent analysis applies to these systems as well. 

The outline of our paper is as follows: We start in section 1 with a general remark 
on the tails of the first passage time distribution in case with independent randomly 
moving searchers. In section 2 we derive an exact result for the distribution P{uj) for 
\E'(r) in Eq. (jl]) with arbitrary /i and for arbitrary A^. In section 3 we discuss several 
exactly solvable particular cases. Section 4 is devoted to the asymptotic analysis of 
P{uj) for large and arbitrary /i. In section 5 we evaluate P{uj) for the exponentially- 
truncated distribution in Eq. (|6]) in case of two searchers. We also furnish here exact 
calculations of P{uj) for two BMs in finite ID, 2D and 3D spherical domains. Next, 
in section 6, we discuss the form of P{u}) for searchers whose first passage time 
distribution is given by Eq. ([6]). Finally, in section 7 we conclude with a summary of 
our results and some generalizations. 
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1. First passage to the target for independent searchers 

We begin with a somewhat evident but conceptually very important remark on the 
tails of the first passage time distribution in case of N independent, non-communicating 
randomly moving searchers. For the original results and a discussion we address the 
reader to Ref. Some other interesting aspects of this model were discussed in 

Refs. dsmsiiie]. 

Let -P/v(t) denote the probability that up to time moment r neither of searchers 
has visited the target, i.e., the target remained non-detected up to time r. For the 
situation under consideration, clearly, 

Pr,{T)=P,{Tf, (7) 

where Pi (r) is an analogous probability for a single searcher. Supposing that the target 
is found as soon as any of the searchers arrives to its location for the first time, we have 
that for N searchers the first passage density \l/Ar(T) is given by 

^N{r) = (8) 

This yields immediately that \l/iv(r) follows (for the parental distribution in Eq. (j4])), 

as r — )■ oo, 

^N{r) ~ (9) 

A remarkable feature of this simple result is that depending on the number of searchers 
N, the distribution in Eq. may have finite moments even if /i < 1, contrary to the 
parent distribution in Eq. (j4]) which does not have any moment for such values of /x. 
For instance, for < yU < 1 and fiN > 1, Eq. has a finite first moment, i.e., a finite 
mean first passage time, for fiN > 2 it has a finite second moment, and generally, for 
fiN > k it has k first finite moments. This is a crucially important advantage of search 
processes involving searchers, which makes the search process more efficient even for 
non-communicating searchers. 

Next, it might be instructive to consider some exactly solvable case, i.e., BM in 
semi-infinite ID systems, and to calculate the iV-dependence of the mean first passage 
time. In this case, the probability that a single searcher has not visited the location of the 
target (the origin) up to time moment r, starting at distance xq from the target, obeys 
Pi(r) = erf (^Xo/y/ADr^ , where erf(.) is the error function. Hence, the first moment of 
the distribution in Eq. ([H]), i.e., the mean first passage time {tn) for an ensemble of TV 
independent searchers is defined as 



dP^ir) 



oo 



The integral in the latter equation is convergent, as we have already remarked, for 
> 3. One verifies that {tj^j) is a slowly decreasing function of the number of searchers 
A^, and, for A^ ^ 1, we find that 
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Therefore, the more searchers one has, the less the mean first passage time is. Curiously 
enough, this essentially ID result, (tat) in Eq. (fTTl) . coincides exactly with the mean 
residence time which BMs spend simultaneously together in a circular disc of radius 
Xq on a 2D plane within an infinite time interval [1^, or with the mean first exit time 
of one of BMs from such a disc |18] . 

Note finally that the mean first passage time may acquire a much stronger 
dependence on and decrease much faster for communicating searchers which share 
the information on the location of the target [11]. Note, as well, that increasing the 
number of "searchers" allows to decrease substantially the time necessary to reach the 
target site on a DNA for proteins which strongly bind to other nonspecific sites acting 
as deep temporal traps [19]. Some properties of the first passage time distribution for 
the event in which A^ random walks appear for the first time simultaneously at the same 
lattice site have been discussed in Ref. [9] within the context of a survival of an evasive 
prey. 



2. General form of P{u) for heavy-tailed first passage time distributions 

Let ^ exp (—Aw) ^ A > 0, denote the moment generating function of the random variable 
u, Eq. ([1]). For arbitrary \&(t), it can be formally represented as an A^-fold integral: 



f f (n *~ exp (-A-^-^^) ,(12) 



Integrating over dri, we change the integration variable ri — )■ u, to get 

' N 

JJdr„^(r, 



Xuj\ / „— Ao; 



7 fnn /Txn / ^ 

auj 



[e ) = I 







vn=2 



X (r2 + . . . + TTv) * (r^^ (r2 + . . . + rjv) j . (13) 
Using next the following integral representation 

POO 

r$(r)= / dpQ{p)e-P\ (14) 



where the kernel Q{p) is some unknown function defined via the inverse Laplace 
transform of the distribution ^(r), we obtain the following general result for the 
probability density -P(w) in case of A^ identic r-variables: 

Pico) = 1- TdXQ f^x] $ (A) , (15) 



with $(A) = (exp(— Ar)) being the characteristic function of the distribution \E'(r). 
Now, for the distribution in Eq.(jl]) we have 

Q{p) = -Y(i7)-P '^'^-^ (2v/^^ , (16) 
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$(A) = ^X^^/^K, (2V^) , (17) 



and 



where Jui') is the Bessel function. 

Substituting the expressions in Eqs. ( IT6|) and (ITTI) into Eq. (ITS!) . we find that the 
probabihty density P{uj) is given by 

^(-) - wr^^ f i^u\ (18) 



= (19) 



independently of a, with 

The result in Eq. (fT8|) defines an exact distribution for arbitrary /i and A^. In 

several particular cases, the integral in Eq. (fT8|) can be performed in closed form: when 
iV = 2 or iV = 3 and arbitrary /i > 0, or when /i is equal to a half of an odd integer, while 
N is arbitrary. We discuss below some of these well as present an asymptotic 

analysis of P{<^) in the limit ^ 1. 



3. Exactly solvable cases for the heavy-tailed distributions 

3.1. Two non- communicating random searchers 

Consider first the case of just two random searchers, N = 2. For this simple situation, 
one readily finds from Eq. (fT8|) that P{uj) is given explicitly by: 

pM = P^c.'^"^(i-^r\ (20) 

i.e., in this case P{uj) is a beta-distribution. This result has been also obtained within 
a different context in Ref. |26j . 

Notice now that, despite its simplicity, the result in Eq. (120|) contains a surprise: 
it has a completely different shape (modality) depending on whether 0</i<l,/i = l 
or yU > 1 (see Fig. [2]). When < yU < 1, the distribution P{oj) has a characteristic U- 
shape, so that the most probable values of u are and 1. Strikingly, the mean (cu) = 1/2 
corresponds in this case to the least probable value of the distribution. This signifies that 
in this case there is no symmetry between two identical searchers and both arrive to the 
target for the first time at distinctly different times. Note that for /i = 1/2 (two BMs 
in ID), P{uj) = l/7rA/cj(l — oo) and hence, the probability P{uj < ^2) that oo attains any 
value from the interval [0, Q] obeys 

2 / \ 

P{uj <n) = -arcsin (^VQJ , (21) 
which is the continuous arcsine distribution. 
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Figure 2. The distribution P{uj) in Eq. ^ for ^ = 1/2, = 1 and ^ = 3/2. 

When n = 1, the distribution P{uj) in Eq. f l2U]) is uniform, which means that for 
either of the searchers its first passage time, relative to the average over the ensemble 
of two searchers, may take any value with equal probability. 

Finally, for /i > 1 the distribution P{uj) is unimodal and centered at u = 1/2. 
This signifies that in this case two searchers will most likely arrive to the target 
simultaneously. 

Therefore, in situations when < /i < 1 (i.e., when \l/(r) does not have already the 
first moment) the first passage times of two searchers will be most probably distinctly 
different. On contrary, for /i > 1 the first passage times will be most probably the same 
and hence, the algorithm resulting in such values of fi will be robust. 

We end up this subsection by noticing that a similar transition was found in Ref. [50] 
in the related Sinai model with a linear drift of strength fi. There the occupation time 
distribution on the positive axis exhibits a transition at yU = 1 and is also given, in 
certain limiting cases, by a beta-distribution, as in Eq. f l20|) . 




3.2. Three non- communicating random searchers 

For three non-communicating searchers, whose first passage times obey a non-truncated 
distribution in Eq. (j4]), we find 

p. ^_ r(2/i)r(3^) 1^^. ^M-i 

^^^)-2^r3(^)r(2^ + i/2)^ ^^"^^ 

X 2Fi(^2/i,3/i;2/i+^;-i^^, (22) 

where 2F1 is a hypergeometric series. The distribution in Eq. f l22|) for three different 
values of // is depicted in Fig. Ei 

Consider first behavior of the distribution in Eq. ( p2|) in the vicinity of a; = 
and u = 1. Recollecting the definition of the hypergeometric series, one immediately 
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Figure 3. The distribution P{uj) in Eq. ([22]) for ^ = 1/2, = 1 and ^ = 3/2. 

observes from Eq. (122|) that P{uj) ~ C2(l — w)^"^ when a; — )■ 1. A httle bit more 
involved analysis shows that -P(w) ~ Ciu;''"^ when — > 0. This means that similarly 
to the N = 2 case, the distribution P{uj) diverges at both edges when fi < 1, and 
most probable values of u are and 1. Note, however, that Ci > C2 and hence, the 
distribution is skewed to the left favoring small values of u. Therefore, for /i < 1 
the most probable situation is that one of three searchers arrives to the location of 
the target much earlier than two others. Clearly, the mean (u) = 1/3 does not have 
any significance [apart, of course, of the fact that this is just the first moment of the 
distribution in Eq. (l22l)]. 

Further on, the distribution in Eq. (!22|) exhibits a qualitative change of behavior 
for /X > 1. Here P{u}) is always a bell-shaped function of the variable u centered at the 
most probable value u = Um- The only difference between the fi = 1 and fi > 1 cases 
is that for the former P{u!) attains a non-zero values at the edges, P{uj = 1) = 1/3 and 
P{uj = 0) = 2, while in the latter P{uj = 0) = P{uj = 1) = 0. It is important to observe, 
however, that Um is always appreciably less than (u) = 1/3. For example, for /i = 3, 
one has Um ~ 0.2719, for fi = 10 one has Um ~ 0.3102 and etc. In fact, Um -t- 1/3 only 
when yU — )■ 00. 



3.3. N non- communicating random searchers for fi = 1/2 and fi = 3/2 

For /i = 1/2 and arbitrary A^, {i.e., for BMs starting from the same point on a 
semi- infinite line), we find the following simple law: 

A^- 1 1 1 
^""^ ~ n ^/^JTT^) 1 - u + {N - lyu' 
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while for /i = 3/2 and arbitrary we get 



c^(l-a; + (iV-l)2u;)=^/'^r(iV-p) 



1 - + (iV - l)2a;y ^^M Vl (iV- l)2a; 

where Ufc(-) are the Chebyshev's polynomials of the second kind. These distributions, 
for a particular case = 10, are depicted in Fig. (jlj). 




Note that for n = 1/2 and arbitrary A^, the distribution P{u!) has a characteristic 
fZ-shape and diverges at the edges a; = and u = 1. The distribution function 
is strongly skewed towards the left edge; for a; — ?■ the distribution behaves as 
P{uj) ~ (iV - 1) w-i/^ while for w ^ 1 we have P(w) ~ (A^ - 1)"^ (1 - i.e., the 

amplitudes are (A^ — 1)^ times different. This means that for A^ ^ 1 an event in which 
one out of a swarm of A^ independent searchers arrives to the target much earlier than 
the rest of searchers is much more probable than an event in which one of the searchers 
arrives to the target location much later than others. 

For fi = 3/2 the distribution P{u) is unimodal and peaked at some value oj = Um, 
which defines the most probable first passage time of a given searcher relative to the 
ensemble average value. Note, however, that Um is always less than the mean (u) = 1 /N. 
Actually, as one may observe from Fig. HI Um is more than two times less than 1/10. 
Note that in this case the mean first passage time for an individual searcher exists, but 
apparently is not a proper measure of the search process since it is different from the 
most probable value. 
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4. Asymptotic large-A^ behavior of P^u) for heavy-tailed first passage time 
distributions 

We consider next the large- asymptotic behavior of the probability density -P(w) in 
Eq. ( ITSl) . To do this, it might be helpful first to look at the Eq. ( fTSl) from a different 
perspective and to realize that P{uj) can be expressed as the probability density function 
P/v„i(|r|) for the position of a random walker in some Rayleigh's-type random flight 
model with a variable flight length in which - (the number of independent searchers) 
- will play the role of "time". 

Consider a Rayleigh's random flight process (see, e.g., Ref. [5l]) in a. d = 2fi- 
dimensional continuum. If d is thought to be an integer, this would restrict the analysis 
to some particular values of fi. However, as noted in Ref. [52] , formal considerations can 
be used to define an isotropic random flight in a space of arbitrary positive dimension 
d, with d not necessarily restricted to integral values. This relaxes any constraint on 
the value of /i. 

Suppose that flights lengths are independent, identically distributed random 
variables and let the probability density function p{r, r') for a transition from r to 
r' be a function of p = |r — r'| only, which means that the process is translationally 
invariant and p{r, r') = p{p). Choose next 

r(2/i) I 

7rA'r(/i)(l+p2 

i.e., a heavy-tailed (but normalizable for any /i > 0) distribution. 

Using Fourier transform technique (see, e.g., Ref. [51]), one readily finds that the 
probability density function of a random walker, starting at the origin, to be at a distance 
p from the origin after TV — 1 such flights, is given by 

1 f°° 

PN^iip) = , , / rfnn^J^_i(pn) S^-i(n), (26) 



PiP) = ^L^'... . ,,2. . (25) 



with E{u) defined by Eq. ([19]). 

Consequently, we find the following relation between the probability density P{oo) 
in Eq. ([18]) and Pv-i(p) in Eq. ([26]): 



u;(l-a;)P(a;)= ^p2''P^..i(p) 



(27) 



P=V (l-t<j)/aj 

Therefore, the probability density of the random variable u, which describes the 
realization-dependent ratio of the first passage time of a given searcher in an ensemble 
of ones, and of the ensemble averaged first passage time, is proportional to the 
probability density of finding a random walker performing Rayleigh's random flights 
with a broad distribution of flight length, Eq. f[25]) . at distance ^yjY^^ZjJJuj away from 
the origin after — 1 flights. 

When A^ is large, the integral in Eq. ( ITS]) is dominated by the behavior of S(m) in 
the vicinity of m = 0. In turn, the latter depends on the value yu: the cases /i > 1, /i = 1 
and < 1 need to be considered separately. 
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^m = -;^{f-i-l + ifi + l)N-J{fi-l + {fi + l)Nf - 8(/i - 1)N 



4.1. The case /i > 1 

In this case, the leading small-w behavior of the characteristic function follows: 

ln(S(M)) ~ -MV4(/i- 1) . (28) 

Plugging the latter asymptotic form into Eq. ( 1T8|) . and performing integration over du, 
we find that the asymptotic large- behavior of the distribution P{uj) is determined by 

- ( ^^V-P (-(^^ - 1)^1 • (29) 

Therefore, in the limit ^ 1 for /i > 1 the distribution P{uj) is always a bell-shaped 
function of w, which approaches exponentially fast when a; — > and as a power-law 
when a; —7- 1. The maximum of P{uj) is located at 

1 

AN 

u — 1 1 , , 

~ . 30 

This substantiates our previous claims that uj^ is always less than (w) = and 
converges to 1/A^ only when yU — > 00. Note that in this case the first moment of ^E'(r), 
i.e., the mean first passage time, exists but it is not representative of the most probable 
behaviour. 

It is instructive now to reproduce the result in Eq. ( 129|) using a different type of 
argument. Notice that Pn{p) in Eq. ( l26l) in the limit N ^ 00 becomes a Gaussian 
distribution in d = 2/i dimensional space of one rescaled variable p/N^/"^: 

^Mp)~]^-p(-(/^-1)^)- (31) 

Inserting the latter expression into the relation in Eq. fl27|) . we recover the result in 
Eq. 



4.2. The case /i = 1 

In this borderline case one finds 

In (H(n)) ~-(l-27 + 21n2-2 ln(u)) u'^/A , (32) 

where 7 is the Euler constant. 

Note that here one has an additional logarithmic factor ln(u), as compared to the 
leading small-u behavior in Eq. fl25]) . Since logarithm is a slowly varying function, we 
can repeat essentially the same argument: Pn{p) in Eq. (126|) converges, as — > 00, to 
a Gaussian distribution of the scaling variable p/ (N log NY^"^: 

1 f P' \ 



Hence, in virtue of the relation in Eq. ( p7|) . we find 

i ( (1-c^ 



- TTTTTTTT^ ^xp ( -^TZT^. ) • (34) 
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This is again a bell-shaped function of lo with the most probable value 
Um ~ l/(2iVlogiV). 

4.3. The case /i < 1 

For /i < 1 the leading small-w behavior of the characteristic function reads: 

ln(S(M)) (35) 

^ ^ 4^r(i + /i) ^ ^ 

so that for -T- 00, -P/v(p) in Eq. (126|) becomes 

where S{-) is the one-sided, 2/i-dimensional stable law with index a = 2fi and scale 
a = r(l — /_i)/r(l + yu) [53]. In consequence, the distribution P{uj) is given, for N ^ 00, 
by 



P a; ~ — S[x = — T77r\ , a = 2u, a = — . 37 

^ ^ Nu^'+^ V A^i/2mV u ^ r(l + /i)/ ^ ^ 

Asymptotic behavior of S{-) has been discussed in detail in Ref. [53]. When x <^ 1, 
S{-) — > const, and hence, 

P(.) ^ (i^. (38) 

z.e., -P(w) diverges as a; — )■ 1. On the other hand, when x ^ 1, i.e., when a; is sufficiently 
close to 0, S{-) ~ 1/x'^'^ [53], which yields 

P{cu) ~ iVw^-^ (39) 

This means that P{u!) diverges when w — )• and, generally, in this domain < /i < 1 the 
distribution has a characteristic ?7-shaped form strongly skewed towards small values of 
u since the amplitudes differ by a factor A^^. Note that for n = 1/2 the distribution in 
Eq. f l37|) becomes the Cauchy distribution so that the simple form in Eq. f l23|) follows 
immediately. 

5. The distribution P{u}) for exponentially-truncated first passage time 
distributions. Two non-communicating searchers 

We turn next to the analysis of the distribution P{u}) in case when the parent first 
passage time distribution \E'(r), Eq. (E]), possesses moments of arbitrary order. 

For two non-communicating searchers and the exponentially-truncated distribution 
in Eq. ([6]), we get the following result for the distribution of the random variable cu: 



= .r.J r^. rTT^ [ 2J—-^ ) . (40) 



2A'2(2vW^) w(l - w) ^ V V ^'^(l - w 

Two remarks are in order. First, one readily notices that P{uj) vanishes exponentially 
fast when w — )■ or — )■ 1 so that P{u = 0) = P{u = 1) = 0. Second, P{oj) is clearly 
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symmetric under the replacement 00 1 — lo. Since here we deal with a truncated 
distribution \E'(r) which possesses the moments of arbitrary order, our first guess would 
be that P{u!) is always a bell-shaped function with a maximum at u = 1/2. To check 
this guess, we expand P{uj) in the Taylor series around u = 1/2: 

2i^2^(4vW^) 



Pioo) 



1 + 



Inspecting the sign of the coefficient before the quadratic term, i.e., 



9 = l-M-2./f ^'-'-''^^^\ (42) 



we notice that 

• For /i > 1, (yf is always negative for any value of h/a so that here the distribution 
P{uj) is a bell-shaped function with a maximum at w = 1/2. 

• For yU = 1, is negative and approaches from below when 6/a — ?■ 00. It means 
that P{oj) is generally a bell-shaped function with a maximum at a; = 1/2, but it is 
becoming progressively flatter when h/a is increased, so that ultimately P{<^) ~ 1 
apart of very narrow regions at the edges for b/a ^ 1. 

• For < < 1 there always exists a critical value ydfJ') of the parameter y = b/a 
which is defined implicitly as the solution of the equation g = 0, Eq. fH2l) . For 
b/a < ydf^), the distribution P{u}) is unimodal with a maximum at w = 1/2. For 
b/a = ydlJ'), the distribution is nearly uniform except for narrow regions in the 
vicinity of the edges. Finally, which is quite surprising in view of the fact that 
in this case ^'(r) possesses all moments, for b/a > ydfJ^) the distribution P{uj) is 
bimodal with a characteristic M-shaped form, two maxima close to and 1 and 
uj = 1/2 being the least probable value. 

We depict in Fig. [5] three characteristic forms of P{u}) for fi = 1/2 and three different 
values of b/a. 

Therefore, for /i > 1 two non- communicating searchers will most probably reach 
the target simultaneously. For < /i < 1 two distinctly different situations are possible: 
if 6/a is less than some well-defined critical value ydfJ'), then most likely both searchers 
will arrive to the target for the first time together. If, on the contrary, b/a exceeds 
this critical value, the event in which both searchers arrive to the location of the target 
simultaneously will be the least probable one. 

5.1. Two independent BMs in a semi-infinite ID system with a bias 

Recall now that the distribution in Eq. ([6]) with fi = 1/2 is exact for a BM in semi- 
infinite ID systems in presence of a constant bias pointing towards the target. In this 
case one has a = Xq/AD, where xq is the starting point and D - the diffusion coefficient. 
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Figure 5. The distribution P{uj) in Eq. ([Ml for /i = 1/2 and different values of b (a 
is set equal to 1). 



and b = AD/v'^, v being the drift velocity. Hence, 2^fajb = Pe = xo\v\/2D is the Peclet 
number (see, e.g., Ref. [33]). Consequently, we can make a following statement: 

Consider two independent, absolutely identical BMs on a semi-infinite line, starting 
at the same point xq, having the same diffusion coefficient D and experiencing the same 
bias F which points towards the origin so that the drift velocity of both BMs is f < 0. 
Then, an event in which both BMs arrive for the first time to the origin simultaneously 
is 

• the least probable if Pe < Pcc , 

• the most probable if Pe > Pec , 

where Pec is the solution of the transcendental equation 

1 = iPe/f'i (43) 

An approximate solution of the latter equation gives Pcc ~ 0.666.... 

Therefore, the mean first passage time to the target might be an appropriate 
measure of the search efficiency for sufficiently large Peclet numbers, but definitely 
is not the one in case of small Pe. In the latter case the sample-to-sample fluctuations 
are significant and the mean value is not representative of the actual behaviour. 



5.2. Two independent unbiased BMs on a finite interval 

We have already remarked that the expression in Eq. (|6]) is a reasonable approximation 
for the first passage time distribution for random motion in finite systems. For a one- 
dimensional bounded interval of length L and for /i = 1/2, the parameter b ~ L^/D, 
while a ^ xl/D, where Xq is the starting point. Consequently, \/ ajb should be xq/L 
and independent of the particles' diffusion coefficient D. This suggests a somewhat 
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strange result that the very shape (or the modahty) of the distribution P{uj) will 
crucially depend on the starting point Xq. 




CO 



Figure 6. Two BMs on a finite interval [0,L] with a reflecting boundary at a; = L. 
The distribution P{u}) in Eq. (|^ for different values of Xq/L, xq being the starting 
point of both BMs. 



To verify such a prediction, we focus now on a BM on a finite ID interval and 
re-examine this question using an exact form of the normalized first passage time 
distribution \E'(r). For the BM with a starting point at xq on an interval of length 
L with an adsorbing boundary at x = and a reflecting boundary at a; = L, \E'(r) is 
given by 



where 



n=0 ^ 



L2 



7r(n + l/2)xo 



(44) 



(45) 



Consequently, the normalized distribution P{u}) in case of two independent, identical 
BMs has the following form 



4 



Alt='^ (xq/L) Alt'^ jxp/L) 
^t^o (^(^ + 1/2)2 + (1 - a;)(m + 1/2)2)2 • 
Using next an equahty 

1 d 1 1 



(46) 



(47) 



^l^du;\l + ^^\l (c^A2 + (l-a;)A^,)2' 
and the following representation of the characteristic function $(A) of the distribution 
^(r) in Eq. dH]), 



<I>(A) = / dr*(r) exp(-Ar) 
Jo 



First passages for a search by a swarm of independent random searchers 



19 





(48) 



2 d 

IT du 



oo 



sm 



{nim + 1/2) f) cosh (7r(m + l/2)y^ (l - f)) 



(49) 



+ V2 cosh (n{m + 1/2) 



The distribution P{co) tends to zero exponentially, P{co) ~ ex^p(—^TXo/2L^/uJ)/uJ^^'^, 
when Cij — 7- 0, precisely in the same way as -P(w) in Eq. ( HOl) obtained for the exponentially 
truncated first passage time distribution. By symmetry, we expect the same behavior 
when a; — 1. In principle, -P(w) in Eq. ( I49l) can be represented in closed form as a 
complicated combination of elliptic function. However, we prefer to proceed with a 
numerical analysis of the rapidly convergent series in Eq. ( l49l) , in order to understand 
whether P{u}) is always a bell-shaped function of u, or undergoes a transition to an 
M-shaped form at a certain value of Xq/L. 

In Fig. [6] we depict the distribution P{u) in Eq. ( l49l) for different values of the 
ratio Xq/L. One notices that P{uj) has a different shape depending whether Xq/L is 
less or greater than the critical value ~ 0.605.... This allows us to make the following 
statement: 

Consider two independent, identical, unbiased BMs starting at the same point xq 
on a finite interval [0, L\ with a refiecting boundary aX x = L. Then, an event in which 
both BMs arrive simultaneously to the origin is 

• the least probable if xq/ L < 0.605..., 

• the most probable if Xq/L > 0.605.... 

From a common sense point of view such a behavior seems a bit counterintuitive 
- indeed, why should two BMs arrive to the target at progressively distinct times the 
closer they are to its location, and should most probably arrive together when they are 
far from it? On the other hand, such a behavior is quite a natural one: Indeed, Xq 
(a^/^) and L (6^/^) define the effective size of the window in which the decay of the first 
passage time distribution is governed by the intermediate power-law tail. The larger is 
this window, the closer we are to the situation described in Section III. Therefore, the 
origin of such a disproportionate behavior of two identical, independent BMs is precisely 
the same as the one behind the famous arcsine law for the distribution of the fraction of 
time spent by a random walker on a positive half-axis [54]: Once one of the BMs goes 
away from the target, it finds it more difficult to return than to keep on going away. 
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Figure 7. A BM on a finite interval [0,L] (with a reflecting boundary at a; = L) 
starting at x ~ Xq. The first passage time distribution ^(r) in Eq. pi)) for L = 100, 
D = 1/2 and different values of xq/L. 



To substantiate this claim, we plot in Fig. [7] the distribution in Eq. (jH]) for three 
different values of Xq and fixed L. Note that all three curves show an exponential 
behavior for both small and large values of r (which, in fact, is an argument in favor 
of our choice of the exponentially truncated distribution in Eq. ([6])). For large r all 
three curves merge which signifies that at such values of r the characteristic decay time 
is dependent only on L. The lower cut-off is clearly dependent only on the starting 
point Xq. Further on, notice that the closer (for a fixed L) the starting point xq to the 
reflecting boundary is, the narrower is the distribution \E'(r). On contrary, the smaller Xo 
is (the closer to the target), the more pronounced the intermediate power-law behavior 
~ becomes (see the dashed line in Fig. [7]). Actually, the fact that the farther 

away the starting point from the reflecting wall is, the broader is the first passage time 
distribution (the BM simply does not "know" that it is in a finite system up to times 
of order ~ L'^/D) has been already discussed in detail in Ref. [33] . 

Consider next some standard characteristics of the first-passage-time distribution in 
Eq. ( 144|) . such as the mean first passage time, the skewness, the coefficient of variation 
and the kurtosis, to see if any of them reflects the transition observed for P{uj) at 
Xq/L ^ 0.605.... Note that the moments of arbitrary order of the distribution \&(t) in 
Eq. ( l44l) can be straightforwardly calculated in closed form: 



where E2k{-) are the Euler polynomials. Consequently, the mean and the variance are 
given by 




(50) 




(51) 
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Figure 8. The coefficient of variation 71,, the skewness 7a and the kurtosis 7e of the 
distribution \E'(r) in Eq. pi)) versus x^/L. Dotted vertical hne defines the critical value 
Xo/L « 0.605... at which the distribution P{uj) changes the modality. 



and 

-4 



respectively. One can readily check that both characteristic properties are monotonically 
increasing functions of Xq/L and do not show any sign of a particular behavior at 
Xq/L ^ 0.605.... Further on, we define the coefficient of variation 



the skewness 



, kl E,ixo/2L)-3El{xo/2L) 

= 'iw ^ V — — ' ^''^ 

7a = i, (54) 
K2 

where is the third cumulant: 

fc3 = (r^)-3(r)(r^)+2(r)^ (55) 
and the kurtosis (coefficient of excess) 

7e = (56) 

with k^ being the fourth cumulant of the distribution in Eq. (jUJ: 

k, = (r^) - 3(r2)2 _ 4(^)(^3^ ^ I2{r)''{r') - 6(r)' . (57) 

In Fig. [8] we plot the coefficient of variation, the skewness and the kurtosis of the 
distribution in Eq. (jH]). One observes a strong variation of these properties reflecting 
an influence of extreme events for sufficiently small values of xq/L. For xq/L > 0.6 the 
variation of 7^ and 7e becomes rather small. However, neither of these properties shows 
a clear demarkation line between different regimes exhibited by P{uj)- 
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We finally remark that in many practically interesting physical problems the 
starting point Xq is not fixed, but the searcher rather starts from some random location 
which is uniformly distributed on the interval. The distribution Pav{uj) appropriate to 
such a situation is obtained by merely averaging P{uj) in Eq. (l49l) over xq, i.e., 

Pav{co) = jj^ dxoP{cu). (58) 

A quick inspection of -P(w) in Eq. (149|) shows that, due to the orthogonality of 
Alf='\xo/L), the averaged distribution Pav{uj) does not depend on u, so that 

Pavico) = 1. (59) 

This seems to be a general property of Pav{oo), associated with the probability 
conservation and thus, as will be checked for further examples, is independent of the 
dimension of space. 

5.3. Two independent BMs in a disc with a reflecting boundary 

Recall that in an infinite 2D system the first passage time distribution has an algebraic 
(with a logarithmic correction) tail with /i = 0, Eq. (I5j). Consequently, we may expect 
essentially the same behavior as we observed in two previous subsections. 




Figure 9. Two BMs in a (2D or 3D) bounded spherical domain of radius L with a 
reflecting (thick black line) boundary. The target is a circle (in 3D, a sphere) of radius 
r centered at the origin. Two BMs start from the same point (a filled circle) at a 
distance Xq from the center and arrive to the target boundary for the first time at time 
moments ti and T2, respectively. 

Consider an immobile target of radius r fixed at the origin of a disc of radius L 
with a refiecting boundary. Suppose next that a BM starts at some point at distance 
Xq from the origin and hits the target for the first time at time moment r. Then, the 
first passage time distribution \E'(r) in such a situation is given explicitly by 
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First passage time distribution ^'(t) in Eq. for = 1/2, L = 100, r = 5 and 
different values of xo/i. 



rD 



n=0 

where Z is the normalization, 

At'\r,xo,L)- 

Uo{x) -- 



Uo{XnXo)UQ{Xnr) 



XlL^UliX^L) - r^U'Mr 
YQ{Xnr)JQ{x) - jQ{Xnr)Yo{x), 
dUo{Xnx) 



dx 



(60) 

(61) 
(62) 
(63) 



while Xn are the roots of the function 

l^o(A„r)Ji(A„L) - Jo(A„r)Fi(A„L), (64) 

arranged in an ascending order, and Yn{-) are Bessel functions of the second kind. Note 
that A„ depends on L and r. 

In Fig. [in] we depict the distribution in Eq. ( 1601) for fixed L, D and r, and several 
values of xq. Note that, apart of the intermediate behavior, which follows a power-law 
with a logarithmic correction ~ l/rlog^(r) (see the dashed line in Fig. [10]), we have 
here essentially the same trend as the one we observed for a BM on a ID finite interval. 
Namely, the closer is the starting point to the location of the target, the broader is the 
first passage time distribution so that the intermediate power-law behavior gets more 
pronounced. 
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Suppose next that we have two BMs starting at the same point some distance xq 
apart of the origin. Then, from Eq. (160|) . we get the following result for the distribution 

^ A'h~'^\r,xo,L) Al^~'^\r,xo,L) 
which can be conveniently rewritten, using the equality in Eq. (147|) . as 

m=0 m \ / 

where $(A) is the characteristic function of the first passage time distribution in Eq. fl60|) 
defined by [33] : 

m = VV; — ;Vx / ^ v ■ ^^r) 

We depict the result in Eq. (pBj) for fixed L and r, and several values of Xq/L. 
One notices that the shape of the distribution P{(jj) is again different, depending on 
the value of the ratio Xq/ L. For xq/ L < 0.643, the distribution has an M-shaped form 
with a minimum at = 1/2, for Xq/-/^ ~ 0.643 the distribution is nearly uniform except 
for narrow regions near the edges, and lastly, for xq/L > 0.643 it becomes unimodal, 
but still it is very broad and has a pronounced flat region around the maximum. This 
signifies that even in this regime the sample-to-sample fluctuations are significant. 

Consider finally the situation when the starting point of two BM is uniformly 
distributed within the space between two concentric circles of radius r and radius L. One 
readily finds, by averaging the result in Eq. (l66l) . that here as in ID one has Pnyiuj] = 1. 
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Figure 12. A BM in a sphere with a reflecting boundary. First passage time 
distribution *(t) in Eq. ([55]) for L = 100, r = 5, D = 1/2 and different values of 

Xq/L. 



5.4. Two independent BMs in a sphere with a reflecting boundary 

We now turn our attention to a Brownian motion in a 3D spherical domain with a 
reflecting boundary. The target is supposed to be a sphere of radius r which is fixed at 
the origin. 

The distribution of the first passage time r of a BM starting at a distance Xq from 
the origin to the surface of the target is given explicitly by 



D 



n=0 



where 



uo{x) = yo{\nr)jo{x) - jo{\nr)yo{x) , 

the derivative MQ(A„r) = duQ{Xnx)/dx\x=r, and 

G{r,L,Xn) = Ljo{XnL)uo{XnL){\nLjo{Xnr) -yo{Xnr)) - 
jo{Kr)yo{Xnr) L-r 



(68) 

(69) 
(70) 



A. 



r^X^ 



(71) 



The set {A^} are the roots of yo{Xnr)ji{XnL) — io{Xnr)yi{XnL) , arranged in an ascending 
order, while j„ (-) and yn(') are the spherical Bessel functions of the first and of the second 
kind, respectively. 

The distribution of the first passage time in Eq. f l60p is plotted in Fig. [T2] for 
three different values of Xq/L. An intermediate power-law ~ r~^/^ is apparent for 
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Figure 13. Two BMs in a sphere with a reflecting boundary. The distribution P{uj) 
in Eq. ([72]) for L = 100, r = 5 and different values of xq/L. 



Xo/L = 0.055, persists for about a decade for xq/L = 0.3 and is entirely absent for 
Xq/L = 0.9, i.e., when the BM starts close to the reflecting boundary. 

From Eq. ( l68l) . we obtain the normalized distribution P{u}) for two independent 
BMs starting at a distance Xq from the origin: 

c,/ N _ 1 A'n~'^\r,xo,L)A'hi~'^\r,xo,L) _ 



n,m=0 \ n ' \ ) m) 

m=0 m \ / 

where $(A) is the characteristic function of ^I/(t) in Eq. (168|) . defined by [33] 



sinh ( i/-^ (L - Xo) ) - \ j^L cosh ( a/^ (L - xq) 
<I>(A) = ^ . (73) 



^0 sinh (^y^(L -r)) - y^L cosh (^y^(L 



This distribution is plotted in Fig. ( !T3|) for fixed L and r, and several values of 
Xq/L. Therefore, also in 3D we find that -P(w) has a different modality depending on 
the value of the ratio Xq/L. For Xq/L < 0.68.., the distribution has an M-shaped form 
with uj = 1/2 being the least probable value. For xq/L > 0.68.., the distribution has, in 
principle, a maximum at w = 1/2 but this maximum is almost invisible so that visually 
the distribution looks more like a uniform one, as compared to the ID case in which the 
maximum is more apparent. 

In case when the starting point xq is uniformly distributed between two concentric 
spheres of radii r and L, we again find Pav{^) = 1- 




Figure 14. A BM starting at a distance xq from the origin of a sphere with an 
adsorbing boundary. The first passage time distribution ^'(r) in Eq. ([71)1 for L ~ 100, 
D = 1/2 and different values of xq/L. 



5.5. Two independent BMs in a sphere with an adsorbing boundary 

Consider finally a geometrically different situation in which two BMs start from the 
same point within a three-dimensional sphere at a fixed distance Xq from the origin but 
now the target is the surface of the sphere. In this case, the normalized distribution of 
the time of the first passage of a BM with diffusion coefficient D to any point on the 
surface of the sphere of radius L from a point at distance xq from the origin is given by 
the series (see, e.g., [61 155]): 

, , 27rD-^, ( Xn\ ( 7r'^n'^DT\ , , 

*(^) = ^ Hi) ( — ) ' (74) 

whose moments of arbitrary order m are defined as 



^ = r(m + 3/2)^UJ ""^"^^l^^J' ^^'^ 
where Bm{-) are the Bernoulli polynomials. The distribution in Eq. f l74p is depicted 
in Fig. [H] for fixed L and D, and different values of the ratio xq/L. Note that here 
the situation is inverse to the one in which the target is situated in the origin - the 
most pronounced intermediate time power- law behavior is observed for xq/ L 1 and is 
absent for small values of Xq/L. The normalized distribution P{uj) in this case has the 
form: 

p, X ^ 4L^ y> / pn+m'^"^ SUl {tIUXq / L) SUl {'KmXQ / L) 




Figure 15. Two BMs starting at a distance a;o from the origin of a sphere with an 
adsorbing boundary. The distribution P{uj) in Eq. (|76l) for L ~ 100 and different 
values of xq/L. 



_ zL2 1 d 

2tixI c^3/2 (1 _ ^)i/2 dy 

where ^3 is the Jacobi theta-function: 

00 

(^3iv,q)= q"^^ exp {2Trimv) . 

m=— 00 

The distribution in Eq. (1761) is depicted in Fig. [T5] for fixed L and different values of 
Xq/L. We observe here a transition from a bell-shaped form with a maximum at u = 1/2 
and an M-shaped form with a minimum at u; = 1 /2 and maxima close to and 1 . The 
transition takes place at xq/L = 0.643.... Averaging Eq. (176|) over the starting point xq, 
we again find that Pav{^^) = 1- 



^3 ^2o(i+M^e--^ 



(76) 



(77) 



6. The distribution P{uj) for exponentially-truncated first passage time 
distributions. Three and more non-communicating searchers 

We turn finally to the situation with more than two searchers whose first passage time 
distribution obeys an exponentially-truncated form in Eq. (Q. For = 3 we find the 
following general result 
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A straightforward analysis shows that -P(w) in Eq. f l78p is always a bell-shaped function 
P((0) 




0.2 0.4 0.6 0.8 1.0 

Figure 16. Three non-communicating Brownian searcliers. -P(tL)) in Eq. (|79p for 
different values of the parameter b (a is set equal to 1). 



for /i > 1. The most probable value Um is, however, always substantially less than 1/3, 
approaching this value only when — oo or 6 — 0. 

The case < yU < 1 is different. Focusing on n = 1/2 {e.g., biased BMs on a 
semi-infinite line), for which Eq. flTSjl simplifies, 

= ''"^l- ^K, I 2J^^-^ 1 , (79) 




we discuss a sequence of different regimes which may be observed when 6/a is gradually 
varied, see Fig. Unj For 6/a <^ 1, P{oj) is peaked at Um ~ 1/3. For larger b/a, Um moves 
towards the origin and P{u!m) decreases (see the inset of Fig. [T7j). For yet larger b/a, Um 
keeps moving towards the origin but now P{um) passes through a minimum and then 
starts to grow. At 6/a ~ 140 a second extremum emerges at a; ~ 0.84 which then splits 
into a minimum and a maximum (see Fig. [T7|) so that P{uj) becomes bimodal. For still 
larger b/a, the minimum moves towards u = 1/2, while the second maximum moves to 
w = 1. 

In fact, the result in Eq. ( !79|) can be straightforwardly generalized for arbitrary N, 
yielding 



The distribution function in Eq. flHOj) shows essentially the same behavior as the one in 
Eq. fl7^ : the only difference is that the critical values of the parameter b/a at which 
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Figure 17. Three non-communicating Brownian searchers. The loci of the extrema 
of P{u!) in Eq. ([7^ . Blue (green) circles define the position of the first (second) 
maximum. Red circles define the position of the minimum of P{io). The inset shows 
the maximal value P(cj,„) vs b (a is set equal to 1). 



Pmax attains a minimal value or when the second maximum emerges depend on the 
number of searchers A^. 

7. Conclusions 

To conclude, in this paper we have studied the distribution P{uj) of the random variable 
CO ~ Ti/r, r = -^""^^fcLi^fc, where r^'s are the first passage times to an immobile 
target by independent searchers, which start their random motion simultaneously 
from the same point in space. Hence, r^'s are independent, identically distributed [with 
distribution \l/(r)] random variables. Since u equals, by definition, the realization- 
dependent first passage time of a given searcher relative to the realization-dependent 
ensemble-average first passage time of searchers, the distribution P{u}) can be viewed 
as a measure of the robustness of a given search algorithm and of the underlying random 
motion (space exploration) , which also probes the validity of the mean first passage time 
of a single searcher as a proper and/or informative measure of the search efficiency. 

We have considered two general forms of \E'(r): The one in Eq. (jl]), which is 
appropriate for search in unbounded domains and is characterized by a power-law long- 
time tail ~ r~^~^, where the exponent /i encodes the specific details of the searchers' 
random motion, and an exponentially tempered form in Eq. ([6]), which is a plausible 
approximation for random motion in finite domains or search assisted by a constant bias 
("smell") towards the target. 

We have shown that for a non-truncated distribution \l/(r) in Eq. (j4]) with /i < 1, 
the distribution of the random variable u has a characteristic tZ-shaped form so that 
the most probable values of u are and 1. For N = 2 the distribution is symmetric 
around u = 1/2 with u = 1/2 being the minimum of the distribution. This signifies that 



First passages for a search by a swarm of independent random searchers 



31 



the "symmetry" between two identical searchers is broken. For N > 2 the distribution 
P{uj) is skewed (by a factor A^^) towards the small values of u and can be expressed, 
in an explicit form, as a one-sided a-stable distribution with a = 2/i. For /i > 1, the 
distribution P{u) has a bell-shaped form but the most probable value of u is always less 
than the mean value (u) = 1/N. The most probable and the mean values of u coincide 
only when /i — > oo. 

For the exponentially truncated "^{t) in Eq. ([6]) the distribution P{u}) always has 
a bell-shaped form for /i > 1. For /i < 1, however, the situation is more complicated 
and interesting. We realized first that for = 2 there exists some critical value of the 
parameter ?/c(yu) = b/a, so that for b/a < ydfJ') the distribution P{uj) has a bell-shaped 
form, but for b/a > ydfJ') it attains an M-shaped form with u = 1/2 being the least 
probable value and two maxima close to the edges of the interval. This signifies that, 
despite the fact that ^'(r) has moments of arbitrary order, two identical searchers will 
arrive for the first time to the target at distinctly different times. Consequently, in 
such a situation the mean first passage is not a proper measure of the search efficiency. 
We remind that such a form of "^{t) is exact for a BM taking place on a semi- infinite 
one- dimensional line with a constant bias directed towards the target. For this physical 
situation, our result implies that two identical BMs, starting at some point Xq and 
having the same drift velocity v < and the same diffusion coefficient D, will most 
likely arrive together to the target (the origin) if the Peclet number Pe = xo\v\/2D 
exceeds some critical value Pcc ~ 0.666. On contrary, if Pe < Pec, an event that these 
two BMs arrive simultaneously to the target location is the least probable event. 

Turning next to an unbiased BM on a finite interval [0, L] we recall that b should be 
proportional to L'^/D and a ~ x^/D. Hence, \/a/b ~ xq/L. This suggests a somewhat 
strange result that for two unbiased identical BMs on a finite ID interval the modality 
of the distribution should depend on how far is the starting point from the reflecting 
boundary. Since "^{t) in Eq. which we have used for the derivation of this result, 
is an approximate form of the first passage time distribution, we have revisited this 
problem using exact forms of ^I/(t) for a BM in ID, 2D and 3D spherical domains with 
a reflecting boundary. We have shown that indeed, the very shape (modality) of the 
distribution P{uj) depends on the ratio Xq/L. 

We have realized that for xq/L < Xcid), where Xcid) ~ 0.61,0.64 and 0.68 for ID, 
2D and 3D, respectively, P{u}) is an M-shaped function of u with a minimum at a; = 1/2 
so that here two unbiased identical BMs will most probably arrive to the target location 
for the first time at distinctly different times. For xq/L > Xc{d)i the distribution has a 
maximum at a; = 1/2 so that, mathematically, the most probable event is that two BMs 
arrive for the first time to the location of the target simultaneously. Note, however, 
that P{oj) is a "bell-shaped" function only in ID (although is still rather broad) but in 
2D and 3D, P{u:) is nearly fiat in an extended region around the maximum and rather 
abruptly vanishes in the vicinity of the edges of the interval. This signifies that here 
sample-to-sample fluctuations are very significant. In case when the starting point xq of 
two BMs is uniformly distributed within the domain (outside the target), we found that 
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P{uj) = 1. This allows us to conclude that in neither of these well-studied situations 
(apart of, with some reservations, ID case with Xq/L ~ 1) the mean first passage time 
of an individual searcher can be considered as a robust measure of the search process 
efficiency. 

We argue that a similar behavior will take place in finite ID systems for fractional 
BM with arbitrary Hurst index H or for a-stable Levy fiights with < a < 1, and, 
more generally, for finite systems of (not necessarily integer) dimension df with fractal 
dimension d^ of random motion trajectories given that a) df < d^ (compact exploration, 
fi = 1 — df/dw < 1) orb) df > d^ (non-compact exploration, fi = df/d^ — I) but 
df < 2duj. 

We have evaluated P{uj) in the case of three and more searchers with the 
exponentially-truncated first passage time distribution in Eq. (j6]). We have shown that 
for /X > 1 the distribution P{u)) is a bell-shaped function of oj for any value of b/a. Next, 
we have demonstrated that for very small values of the ratio b/a the distribution P{oj) 
has a bell-shaped form with a maximum close to 1/N. Further on, we have predicted the 
following sequence of regimes which can be observed upon gradually increasing b/a (the 
starting point Xq of searchers is moving towards the location of the target): the most 
probable value Um of u moves towards the origin. The value of the maximum, P{ujm) 
first decreases, passes through a minimal value and then starts to increase. At a certain 
threshold value of b/a, a second extremum emerges in the vicinity of u = 1, which then 
splits into a minimum and a maximum so that the distribution P{uj) becomes a skewed 
M-shaped one. This signifies, as well, that for such a situation the mean first passage 
time of a given searcher is not a representative characteristic of the search process. 

As a final observation, we note that one may encounter a power-law distribution 
with a more abrupt truncation, compared to the exponential function in Eq. ([6]), say, a 
bounded power-law (see [31]) or a power-law tempered from both sides by a Gaussian 
function: 



One finds that P{uj) in Eq. (!82|) is always a bell-shaped function for f-i > 1. For /i < 1, 
depending on the value of z = 2a/b, it may have a unimodal, or a three- modal form. 
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(81) 



In this case, the distribution of the random variable u is given by 




(82) 
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Figure 18. P{(^) in Eq. for a faster than exponential truncation of the 

intermediate power-law behavior. Here the parameter z ~ 2a /b. 
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